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STABILITY OF STREAMWISE VORTICES 


By 

M. R. Khorrami, 1 C. E. Grosch, 2 and R. L.-Ash 3 

ABSTRACT 

A brief overview of some theoretical and computational studies of the 
stability of streamwise vortices is given. The local induction model and 
classical hydrodynamic vortex stability theories are discussed in some 
detail. The importance of the three-dimensionality of the mean velocity 
profile to the results of stability calculations is discussed briefly. 

In this study the mean velocity profile is provided by employing the 
similarity solution of Donaldson and Sullivan. The global method of Bridges 
and Morris was chosen for the spatial stability calculations for the non- 
linear eigenvalue problem. In order to test the numerical method, a second 
order accurate central difference scheme was used to obtain the coefficient 
matrices. It was shown that a second order finite difference method lacks 
the required accuracy for global eigenvalue calculations. Finally the 
problem was formulated using spectral methods and a truncated Chebyshev 
series. 
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1. INTRODUCTION 

The study of the stability of streamwise vortices is an important and 
challenging problem. The breakdown of leading edge vortices on delta wings, 
which severely reduces lift, and the very stable wing tip vortices shed from 
large commercial aircraft, which determine the flight frequency at airports, 
are two classical examples. Meaningful experimental and computational work 
is rare due to the complex nature of the phenomena, and theoretical studies 
are incomplete. The complexity stems from the three-dimensional nature of 
the problem and hence the necessity of finding three-dimensional mean veloc- 
ity profiles which are solutions of the Navier-Stokes equations. Most of 
the stability calculations to date have used Lamb's vortex [Lamb, 1945J, and 
in some cases a Long's vortex [Long, 1958] or solid body rotation super- 
imposed on Poiseuille flow in a pipe as the basic unperturbed flow. Here we 
give a brief overview of some of the more important theoretical and computa- 
tional studies of these flows. This review is by no means exhaustive, but 
the cited papers are those which we have judged most relevant to the present 
research. 

The theory of vortex stability or vortex deformation has progressed 
in two distinct directions. In the first case one uses a localized-induc- 
tion model for an inviscid incompressible vortex [Hama, 1962]. The second 
approach uses the classical hydrodynamic stability theory for rotating 
fluids and can be attributed to a paper of Howard and Gupta [1962] which 
contains a generalization of a stability criteria put forward by Lord 
Rayleigh [1916]. In this summary, we will discuss the development of both 
theories starting with the localized-induction model. 

The deformation of a curved vortex filament under its own influence was 
first calculated by Hama [1962] using a localized-induction model with the 
assumptions that the local radius of curvature is much greater than the 



core radius and that the vortex filament is unaffected by the events at 
infinity. He used several initial two-dimensional shapes and found that in 
every case the region near the vertex of the curve would lift off of the 
plane of the initial curve causing a three-dimensional helical deformation 
to take place. The helical wave, which rotates opposite to the circulatory 
motion of the vortex filament propagates away from the vertex and its 
amplitude increases as the vortex lifts. 

Betchov [1965] extended Hama's work and derived two equations describ- 
ing the curvature and torsional effect of the vortex filament. After lin- 
earizing the equations, he found that the vortex became unstable if the 
perturbation wavelength was > 2irR, with R the local radius of curvature. 

This type of perturbation is generally referred to as a long wavelength 
instabi 1 ity. 

The work of Betchov was extended by Widnall [1972] who considered the 
effect of sinusoidal disturbances on a helical vortex filament of finite 
core radius and infinite extent. She found that the helical vortex filament 
had three distinct modes of instability: a very short wave instability mode 

where the local radius of curvature of the helix was the characteristic 
length; a low wave number mode which was found to be due to the local-induc- 
tion (this instability is the same mode that was identified by Betcnov); and 
the mutual inductance mode. The latter mode sets in whenever the successive 
turns of the helix get very close to each other. Finally, Widnall found 
that increasing the vortex core size: (1) reduced the amplification rate of 
the large wave instability; (2) increased the amplification rate of the 
mutual-inductance instability; and (3) decreased the wave number of the 
short-wave mode. 
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The hydrodynamic stability of a vortex was first studied by Lord 
Rayleigh [1916]. In his classical paper he was able to show that a neces- 
sary and sufficient condition for stability of an inviscid vortex without 
axial velocity was that the square of circulation must increase outwards. 
Much later, Howard and Gupta [1962] were able to give a sufficient condition 
for the stability of an inviscid vortex with axial velocity (subject to 
axi symmetric disturbances). They found that for stability the local value 


J(r) 


1 

r 3 



> 


1 

4 


everywhere. Here r is the circulation and W is the axial velocity. 

It is obvious from the above criterion that axial shear has a destabi- 
lizing effect. For non- axi symmetric disturbances, Howard and Gupta found 
that a sufficient condition for stability was 

„*♦ - V * - i [a ^ + n d(v/r) ] 2 > 0 
r 2 dr 4 dr dr 

where 

, _ 1 d(r 2 V 2 ) 

<P = —r 

r^ dr 


and V is the tangential velocity, and a and n are the axial and azimuthal 
wave numbers, respectively. They stated that the above criterion is always 
violated for sufficiently small a, and while this does not imply instabil- 
ity, it has led them to suggest that no general necessary and sufficient 
criterion is obtainable. 

Using the equation for the radial amplitude of disturbances derived by 
Howard and Gupta, Pedley [1968] has shown that for a very small Rossby 


3 


number (e = — 2_ , where n is the angular velocity) the flow is unstable to 
2ftro 

n on- ax i symmetric inviscid disturbances of sufficiently large axial wave- 
length. He found that, although both solid body rotation and Poiseuille 
flow in a pipe are stable with respect to infinitesimal disturbances, the 
superimposed combined flow (which is a helical vortex) is highly unstable 
witn the negative azimuthal wavenumbers as the dominant unstable modes. 

These disturbances are helical in shape, wrap around the helical vortex in 
the opposite direction and are either standing waves or are traveling up- 
stream. The growth rate of the most rapidly growing disturbance is 2efi, 
where e is the Rossby number and ft is the angular velocity of the pipe wall. 
In a follow-up paper, Pedley [1969] has shown that for viscous, rotating 
Poiseuille flow the critical Reynolds nimber. Re , has a value of 82.9 

v 

corresponding to n = 1. The critical Reynolds number increases as the value 
of azimuthal wave number, n, increases. The disturbances are stationary 
relative to the rotating frame of reference, and as the Reynolds number 
increases, the wave nunber of the most rapidly growing disturbance also 
increases. 

Maslowe [1974] has studied the same flow field as that of Pedley [1968] 
without making an assumption as to the magnitude of the Rossby number. He 
has shown that the most unstable modes have negative azimuthal wave numbers. 
They spiral in the same direction as the basic flow rotation (note this is 
in contradiction to Pedley's finding and is caused by different interpreta- 
tions of what direction a wave with negative azimuthal wave number spirals) 
but propagate upstream in the axial direction with an axial phase speed 
0(e -1 )* Also, the amplification factor and the axial wave number of the 
fastest growing disturbance peaked at a finite Rossby number e of 0(1) while 
the growth rate showed little variation with n for |n|>2 at finite values of 
e . 
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The stability of the mean velocity profile of a trailing line vortex 
(Batchelor [1964]) was studied by Lessen, Singh, and Paillet [1974] in a 
study of the inviscid stability of swirling flows with respect to infinites- 
imal non-axi symmetric disturbances. It was found that negative azimuthal 
wave numbers are destabilized by the addition of swirl. Lessen, et al. 
discovered that the stability of the vortex is very dependent on the value 

P 

of q, where q is given by __2 , with U and r the scaling factors for the 

Ur s s 

s s 

velocity and radial coordinate, respectively and r 0 the constant circulation 

at large radial distance ro . All wavelengths appear to become damped, and 
the flow completely stabilized at a value of q slightly greater than 1.5. 
Using a finite-difference method. Duck and Foster [1980] have solved exactly 
the same proDlem as that of Lessen et al. [1974] and obtained similar 
results. However, they found that the number of unstable modes is directly 
dependent on the number of grid points, N, so that as N increases the numoer 
of unstable modes increases as well. One has to be skeptical about this 
finding since the number of unstable modes should not depend on the grid and 
grow without bound. 

Following their earlier work. Lessen and Paillet [1974] performed a 
viscous stability calculation for a trailing line vortex. They obtained 
similar results to those of Pedley [1969] which indicated that the critical 
Reynolds number increases as |n| increases. Their calculations of neutral 
stability curves showed that the values of the critical wavelength and 
critical Reynolds number are not very sensitive to the exact value of q; 
and at large q all of the unstable modes are stabilized for any wavelength 
and Reynolds number. 
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Although these calculations are valuable, they do not constitute a 
systematic study of the effects of different boundary conditions and mean 
velocity profiles on the stability of the vortex. The need for such a sys- 
tematic study is the rationale behind the present research. In this study 

the effect of such parameters as mean velocity profile, vortex Reynolds 
r dP 

number — , and axial pressure gradient — on vortex stability will be con- 
v dz 

sidered. We are going to examine the spatial stability of a rather general 
class of laminar incompressible vortices in the context of linear tneory. 

2. MATHEMATICAL FORMULATION OF THE STABILITY PROBLEM 
2.1 Mean Velocity Profile 

The similarity solution for porous pipe flow due to Donaldson and 
Sullivan [1960] has been selected as the mean velocity profile. Although 
one might argue that this similarity solution is for confined flow and its 
relation to the unconfined flow of longitudinal vortices is not an obvious 
one, the wealth of different solutions possible (ranging from a single cell 
vortex to multiple cell vortices) make the selection a reasonable choice. 
This multiplicity in the core of some vortices has been shown experimental- 
ly, (Adams and Gilmore [1972]). These experiments also indicate that in 
order for a vortex to become unstable, large gradients in the axial direc- 
tion must develop (Graham and Newman [1974], Leuchter and Solignac [1983]). 
They have shown that the transition of the axial velocity profile from a 
strong jet to a strong wake plays a dominant role in the loss of stability 
of these vortices. Therefore any velocity profile selected for study should 
show similar behavior. 

The mean velocity profiles which we have selected for study are of the 

form 
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U = U(r) 

V = V(r) 

W = zW(r) 


where, U, V, and W are radial, tangential, and axial velocities, respective- 
ly, and r is the distance in the radial direction. Note that the flow 
undergoes linear acceleration in the streamwise direction. This interesting 
feature, we believe, may strongly effect the spatial stability of vortices. 

Figures 1 thru 3 show three sample mean velocity profiles obtained 
using Donaldson and Sullivan's similarity solution. The similarity of the 
tangential velocity in Figs, lc and 2c to that of Lamb's vortex is note- 
worthy. Figures 2b and 3b represent the mean velocity profiles of two 
multiple cell vortices with jet like and wake like axial velocity, respec- 
tively. It is important to note that the transition of a vortex core from a 
jet flow to a wake type can be simulated step by step by changing a few 
parameters in the similarity solution. 

2.2 Perturbations of the Governing Equations 

Cylindrical coordinates (r, 0, z) have been chosen as the coordinate 
system. The equations of motion in cylindrical coordinates are: 

continuity 
r-momentum 


1 3 , , v . 1 3v' 3w' n 

_ (ru‘ ) + _ + = 0 


U) 


r 3r 


r 30 


3Z 


3U ' , , 3U* V' 3u' t . 3u* V' 2 

+ u + + W - 

3t 3r r 30 3z r 
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- 1 9 P' * (r,2 . U ' 2 3V' 

- - — + v (v z u' - — - J 


p 9r 


r 2 r 2 30 


0 -momentum 


£v| 3v' V' 3 V 1 A , 3v ' llV 

+ u + + w + 

3t 3r r 39 3z r 


1 3 p ' , . v' 2 3 u'-> 

= - — + v ( V 2 v‘ - _ + _ ) 


pr 30 


r* r 


2 30 


z- moment urn 


where 


3w^ +u , + «, Sw^ . _ 1 9p^ +v 72w , 

St 


>' 3w' 

3r r 30 


3z 


P 3Z 


, !i_ + i I_ + L + 

3r 2 r 3r r 2 30 2 3z 2 


( 2 ) 


(3) 


(4) 


The perturbed flow variables in these coordinates are assumed to be of 
the form: 


u' = U(r) + u 

v' = V(r) + v 

w' = W(r,z) + w (5) 

P* = P(r,z) + p 

where u, v, w, and p are the perturbations in the velocities and the pres- 
sure. These quantities are assumed to have a helical wave form: 

fu, v, w, p} = fiF(r), G(r), H(r), P(r)> e l(ctZ + n0 " ut) (6) 
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Here, F, G, H, and P are the perturbations amplitudes which are func- 
tions of radial position only, a is the wave number in the axial direction, 
n is the wave ntmber in the azimuthal direction, and w is the frequency of 
oscillation. For a single valued solution, n must be an integer. The case 
of n equal to zero corresponds to that of axisynmetric disturbances and 
positive and negative integers represent different directions of propaga- 
tion. For positive n, the helical wave spirals in the same direction as the 
circulatory motion of the vortex while a negative value indicates that the 
wave spirals opposite to the circulatory motion. 

Substituting into the equations of motion and neglecting the second 
order terms, we obtain the small disturbance equations 


continuity 


3u 

u . 

i 

3v . 

. 3w 

+ _ + 


— T = 

3r 

r 

r 

90 

3Z 


(7) 


r-momentum 

ll + U 11+ u — + 1 3u + u 3u . 2Vv , , jig 
it Jr dr r 90 3z r p 9r 


+ 




r 9r r 2 


3 2 u 

30 2 



u _ 2 3v j 
r 2 r 2 30 


(8) 


0 -moment urn 


3v „ 3 v , 

— + U — + u 

3t 3r 



V 

r 


3v + w 3v + Vu + _vU 
30 3 z r r 


1_ ip 

pr 30 


+ v [ 


3 2 V 




3r 2 r 3r 


3 2 v 

30 2 


3 2 V 
3Z 2 


+ Lli] 

r 2 30 


(9) 
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z-momentum 


3w 3W 3w V 3w „ 3w 3W 1 3p 

— + u — + U — + — — + W — + w — = - ti 

at ar ar r 88 8z az p az 

+ v [ 32w + * aw + * a 2 w + a 2 wj 

3r 2 r 3r r 2 30 2 3z 2 


( 10 ) 


The boundary conditions at the outer wall are 


at r = R 0 

and on the centerline 1 


at r = 0 


u = 

0 

V = 

0 

w = 

0 

p = 

0 


lim 

O 

II 

♦;rl 

r-*-0 

39 

lim 

1- 

II 

o 

r-*-0 

30 


(ID 


(12) 


where q is the total velocity vector. 

Substitution of the perturbation Eq. (6) into the linearized momentum 
equations results in four equations for the disturbance amplitudes. These 
equations in non-dimensional form are: 


continuity 

F ' + I + Hi + a H = 0 (13) 

r r 


x We are very grateful to Dr. M. Y. Hussaini of the Institute for Computer 
Applications in Science and Engineering for mentioning this simple yet very 
powerful way of determining the boundary conditions on the centerline. 
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r-momentum 


ill + i [U - —]r + [« + i - 01 - ctW + i- (Hi 
Re Re r dr r Re r 2 


* * V )] F * [^- 

r 2 Re r 2 r 


£1] G + P' = 0 


( 14 ) 


e -momentum 


z-momentum 


i_ G" + [ U - - 1 ] G' + [-iai + 1H1 + iaW + .1 + -L (— 

Re Re r r r Re r 2 


+ a 2 + — ) ] G + f 1 — + — + — ] F + * 0 (15) 

r 2 dr Re r 2 r r 


— H" + [U - J_ ] H‘ + [-i<u + iaW + ml + — + 1_ (Hi 
Re Re r r 3z Re r 2 


+ a 2 )] H + i m F + iap = 0 
3 r 


(16) 


where Re is the Reynolds number based on pipe radius Ro and prime denotes 
differentiation with respect to the radial coordinate. The boundary con- 
ditions at r = R 0 are: 


F (Ro ) = G(Ro) = H(Ro) = 0. 

P'(Ro) = Normal momentum at Ro. 

The derivation of conditions on the centerline is not simple or 
straight forward and deserves a fuller explanation. In expanding the con- 
ditions (12), we need to consider only the perturbation part of the velocity 
since the mean velocity profile is independent of the azimuthal direction. 
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Assuming that the total perturbation velocity field is represented by v, we 
have: 


or 


1 im 
r+0 


3v 

99 


3v 

36 


— (u e r + v e 0 + we z J 


3u > 



+ 



de 



0 de ae 


e + 
z 



de 


But in cylindrical coordinates (Appendix 2, Batchelor [1967]) 



de 


de 


r 


de 


+ 




de 



Substituting for the velocities from (6) and evaluating the deriv- 
atives, we deduce 


1 im 
r+0 


— = (-nF-G) e r + (iF + inG) e ft + inH e = 0 
36 0 L 


In order for the equality to hold, each component of the resultant vector 
must be zero. Similarly, the limiting process applied to the pressure leads 
to 
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In summary, we have at r = 0, 


nF + G = 0 
F + nG = 0 
nH = 0 
nP * 0 


(17) 


The above conditions depend on the value of the azimuthal wave number, n, 
such that 


If n = 0 


If n = ± 1 
If |n| > 1 


F(0) = G(0) 

= 0 

(18) 

H(0) & P(0) 

are finite 

F(0) ± G (0) 

= 0 

(19) 

H(0) = P(0) 

= 0 

F(0) = G (0) 

= 0 

(20) 

H(0) = P(0) 

= 0 . 


Due to the extreme complexity of the Eqs. (13)- (16), no attempt was 
made to reduce the above system into a single, sixth order ordinary differ- 
ential equation. We believe that this would be of no great value in solving 
the problem numerically. 


2.3 Proposed Numerical Scheme 

The system of Eqs. ( 13) - ( 16) can be represented in a compact vector 


form as: 


1(a) £ = 0 (21) 

with 



and 

L(a) = D Q a 2 + 0^ + D 2 (22) 

Do, Di , and D 2 are discretized versions of the differential operator 
matrices. This type of operator (where the eigenvalue appears nonlinearly) 
is called a Lambda-matrix (Lancaster [1966]). 

It was decided to use the method of Bridges and Morris [1984] to solve 
this nonlinear matrix eigenvalue problem. The advantages of choosing their 
method are threefold. (1) It is a global eigenvalue scheme. This is clear- 
ly an advantage because the availability of an initial guess for a shooting 
method is not needed. Of course these global eigenvalues can be refined 
until a desired accuracy is reached, by any local iterative scheme once they 
are obtained. (2) The method is robust. In some cases infinite eigenvalues 
are encountered, and in such instances, the remainder of the eigenvalues are 
obtained without any difficulty. (3) Most of the subroutines needed are 
available in any mathematical software library. In particular this is true 
for the complex version of the QZ routine which is the heart of the present 
scheme. 

The method can be described quite briefly. First the discrete operator 
L(a) is factored such that 
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L(o) = [Doo + D 0 Y + Di] [la - Y] + DoY 2 + DiY + D 2 . 


(23) 


where I is the identity matrix and the Y matrix is referred to as the right 
solvent (Gantmacher, 1959). Now, in order that (23) be consistent with 
(22), one must have 


(24) 

(25) 

Multiplying Eq. (25) by a and subtracting from Eq. (22) yields 

aDo Y = D 2 . 

Since a is an arbitrary scalar multiplication factor, in general 


[I a - Y] = 0 
or 

[D oa + D 0 Y + Dj] = 0. 


aDoY * D 2 (26) 

and Eq. (24) must be satisfied. Therefore from the above factorization, it 
is obvious that what we must find is the root matrix of the matrix poly- 
nomial 


DoY 2 + DiY + D 2 = 0 

This can easily be accomplished through an iterative procedure. Then a 
subset of the eigenvalues of 1(a) is given by the set of eigenvalues of Y 
obtained from (24). A flow chart of the numerical scheme is given in Fig. 4 
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which shows the steps needed to find a single solution. For a further 
description of the method and its implementation, the reader is referred to 
the paper of Bridges and Morris [1984]. 

3. TEST PROBLEMS 

3.1 Calculation of the Matrices and Testing of the Algorithm 
Prior to the implementation of the difference scheme, toe algorithm was 
checked using a problem with known eigenvalues as a test case. This problem 
was taken from Gregory and Karney [1969], and along with its eigenvalues is 
given in Table 1. The correct eigenvalues were obtained and are tabulated 
in Table 2. Finding the right values assured us that the coding was free of 
logical error. 

It is well known that spectral methods are efficient methods for solv- 
ing eigenvalue problems. Metcalfe and Orszag [1973] have shown that for 
calculation of the linear stability of these types of swirling flows spec- 
tral calculations are very accurate. Because there are no results available 
with which to compare our results (for the stability of velocity profiles 
due to Donaldson and Sullivan), it seemed desirable to solve the above 
matrix eigenvalue problem using two independent formulations. It was 
decided that in order to obtain an approximate estimate of the eigenvalues, 
finite differencing of the operator matrices would be employed. These 
approximate values could then be used to cross check the eigenvalues found 
by a spectral method. This was done but peripheral problems were encounter- 
ed relating to difficulties with convergence (actually, lack of convergence) 
of the solution. A new test problem was chosen in order to check the accu- 
racy of the finite difference approximation. This problem is that of 
Poiseuille flow in concentric cylinders and is described in section 3.2. 

After careful examination, we found that the lack of convergence of the 
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Table 1. Matrix used as a test problem. 


5 + 9i 5 + 5i 

3 + 3i 6 + lOi 

2 + 2i 3 + 3i 

1 + i 2 + 2i 


n 


-6 - 6i 

-7 

- 7i 

-5 - 5i 

-6 

- 6i 

-1 + 3i 

-5 

- 5i 

-3 - 3i 


4i 


Eigenvalues: 

Xi = 1 + 5i 

A 2 = 2 + 6i 

A3 = 3 + 7i 

A 4 = 4 + 8i 


Table 2. Calculated eigenvalues using the method of Bridges and Morris. 


REAL C 

IMAG C 
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6.000000 

4.000000 

8.000000 

3.000000 

7.000000 

1.000000 
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difference solution was dependent on a number of factors (see Appendix A) 
namely, boundary conditions on the centerline, number of grid points, and 
second order accuracy. Therefore, a model problem had to be chosen so that 
an analysis of some of the above mentioned factors was possible in a short 
time and at the same time the problem had to have been solved oy other work- 
ers to enable a comparison. 

The temporal stability of Poiseuille flow in an annulus was picked to 
be the model problem used for testing. 

3.2 Temporal Stability of Poiseuille Flow in Concentric Cylinders 

For this problem there are no effects of the centerline boundary con- 
ditions because the no-slip condition applies at both boundaries. Further- 
more, the continuity equation along with the pressure terms were staggered 
and evaluated at the set of points { r j + ^/2^ w ^ e the momentum equations 

were evaluated at the set of points {r.}. Here the index j takes on values 

<1 

between zero and N, where N specifies the total number of nodes. The grid 
points were distributed uniformly over the gap distance and no attempt was 
made to concentrate or stretch them. Therefore, the only factors influenc- 
ing accuracy that we had to deal with were the nimber of grid points and the 
second order central differencing scheme. 

3.2a Mean Velocity Profile 

The mean velocity profile for Poiseuille flow in an annulus is 


U = 0 
V = 0 
W = W(r) 
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where, in non-dimensional form, W is given by 


W _ 


W 


M 


1 - r 2 + r^ zn r 2 
M 

1 - r 2 + r 2 u,n r 2 
M M M 


( 27 ) 


The non-dimensional ization has been done with respect to the maximum 
velocity and the half gap distance (see equations 9B and 4B in Appendix B) . 
A detailed derivation of this velocity profile can be found in Appendix 
B. 

3.2b Governing Equations 

It is obvious from the form of the mean velocity profile that the 
equivalent of Eqs. (13) - (16) is much simpler due to the absence of both 
radial and tangential velocity components. The equations are: 
continuity 

— + F' + i!i+aH = 0 (28) 

r r 

r-momentum 

- ]_ F" - J_ F' + [a) - aW + 1_ (Hi + a 2 + 1_)] F + l2n G + P* = 0 (29) 

Re Re r Re r 2 r 2 Re r 2 

0-momentum 

- !_ G" - — — G' + [ -iw + iaW + - — (— + a 2 + — )]G 

Re Re r Re r 2 r 2 

+ 2n F + H P = 0 (30) 

Re r 2 r 
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z-momentum 


- 1_ H" - _L H' + [-io» + ictW + L_ f!ll + a 2 )] H + i F + iaP = 0. (31) 
Re Re r Re r 2 dr 


Employing a second order accurate central differencing technique for the 
derivatives, the discretized equations are 


F . + F . , F. - F. 


) >1 + J ii + 1 (G, + S, ,) *1(H, + H. ,) 

?r . . Ar J ^ 1 9 J J 1 


2r j-l/2 Ar 2r j-i/2 


+ ywPj = 0 


(32) 


i_ [ F j-n ~ 2f j * F j-i ] + -i_ ^>1 ~ F j-ij + j 


Re Ar 2 


Re r . 2Ar 


(D - aW. 
J 


+ — + a 2 )] F. + J il l G. + [lit 1 . i] = 0 (33) 

Re r< 2 J Re r,- 2 J Ar 


* * 
P,.,- P, 


Re r. 2 

J 


Re r- 2 J 

J 


i G . ~ 2G . + G . - , G - G . . 

— r — - - —] + ll + [- ito + iaW. 


Re Ar 2 


Re r . 2Ar 

J 


1 ,n 2 +l „ 2n _ in r * *, 

— C r + ° 2 )] G i + r F i + ^ p i+l + p i] = 0 

Re r. 2 J Re r- 2 J 2r. J 1 J 

J J J 


(34) 


i H. + , - 2H. + H. . . H. +1 - H. . 

_ [_ii 2 2li] + ll — [ -2.1 2ll] + [- ito + i aW. 

^ Drt 0 j 


Re r . 2Ar 

J 


* « 2 )] Hj ♦ i fi, F, + (P j+1+ P.) - 0 

Re p . 2 J rl»* 3 J 9 J + -*- J 


dr 


ia * * 

2 


(35) 


where asterisks denote mid-cell values. The boundary conditions 
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No slip 


No penetration 


G(a) = G 0 = 0 

H(a) = H 0 = 0 

' G(b) = G n = 0 

H(b) =H n -0 

F(a) 5 F 0 = 0 

F(b) = F N = 0 


where a and b are the radius of inner and outer walls, respectively. Note 
★ 

that the term ywPj has been added in Eq. (32). This makes the coefficient 

matrix of omega non-singular. Otherwise it would have N rows with zero 

9P 

entries. This term is obtained by adding y — to the continuity equation, 

at 

which is why it is called an artificial compressibility factor (Malik and 
Poll [1985]). The parameter, y, is chosen to be a very small number and, in 
the present case, a value of 10‘ 12 was used. This term causes large values 
for some of the unimportant eigenvalues; however, its affect on the desired 
eigenvalues is negligible. Experimentation with the value of y with the 
present code indicated that this was indeed the case. 

The above system can be expressed in the following way: 


where 


[K-(uM] [X] = 0 



The eigenvalues are then obtained from the requirement: 
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[M _1 K— to I ] = 0 . 

-1 * 

Here we have made sure that M exists by adding the term ytoPj (as mentioned 
previously) to the continuity equation. The remaining procedure is straight 
forward and simple. Any Gauss-Oordan algorithm with pivoting strategy can 
be employed to invert matrix, M, and then using a standard QR routine, the 
eigenvalues of the matrix M * K can be evaluated. 

3.3 Results 

The above system of equations and boundary conditions were solved on 
the vector processor computer (VPS-32) at NASA Langley Research Center for 
the narrow gap case. It is well known that as the gap distance shrinks, for 
fixed inner radius, the eigenvalues approach those of plane Poiseuille flow. 
Therefore, the narrow gap calculation was performed so that the results 
could be compared to the accurate eigenvalues of pi ane-Poiseui 1 le flow 
reported by Orszag [1971]. All of the calculations reported here correspond 
to the axi symmetric type of disturbances where the azimuthal wavenumber, n, 
is zero, and the gap distance is 0.01. The results of these calculations, 
along with Orszag' s results, are tabulated in Table 3. 

Initially, the calculations were started with 20 grid points and no 
satisfactory results were obtained. That suggested that the number of grid 
points had to be increased. After increasing the grid points, many of the 
eigenvalues found by Orszag were obtained, but with only modest accuracy. 
The number of nodes was increased further by fifteen percent from 140 to 
160. However, improvement of the accuracy of the eigenvalues was slow with 
the side effect of almost doubling the cost of each computer run. 

We concluded that a second order accurate finite difference scheme is 
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10000. a = 1, Re = 10000. 



only accurate enough for obtaining rough estimates for this type of global 
eigenvalue calculation and, therefore, the method was abandoned. Hence, a 
spectral method has been selected. 

4. SPECTRAL METHOO 

A Chebyshev collocation approach was chosen. Truncated series of 
Chebyshev polynomials were employed in expanding the flow variables. 

Because the problem is non-periodic, this form of expansion has the advan- 
tage of eliminating terminal discontinuities. In addition, the expansion 
exhibits a rapid convergence rate as the number of terms increases and clus- 
ters the collocation points (see Eq. 43) near the boundaries. 2 The 
Chebyshev polynomials are defined on an interval of (-1, 1) by 

T k (£) = cos [k cos- 1 ?] (36) 

Because the physical range in this problem is (0, 1), a simple transforma- 
tion is made from the physical variable r to a transform variable £ by 

? = 1 - 2 r (37) 

where 

-1 < ? < 1 . 

The mean velocity profile must be transformed similarly in order to be com- 
patible. A sample case is shown in Figs. 5a-5c which correspond to Figs, 
la-lc. 

2 For a thorough exposition of Chebyshev approach, the reader is urged to 
consult the excellent treatises of Gottlieb and Orszag [1977], Gottlieb, 
Hussaini and Orszag [1983], and Hussaini, Streett and Zang [1984]. 
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Let 


then 


C = cos 6, 

T k U) = cos ke. (38) 


The periodic nature of Chebyshev polynomials is apparent from the above 
identity. The governing equations are then satisfied exactly at collocation 
points which are the roots of the Nth degree Chebyshev polynomial T^U)- It 
must be mentioned that it is only at this set of points that one obtains 
spectral accuracy (Gottlieb and Orszag [1977]). These points are defined 
by 


where 


TTJ 

r . = cos — 
J N 


(39) 


j • 0,1,2 N 


with j = 0 and j = N corresponding to the centerline and wall boundary con- 
ditions, respectively. 

An interpolant polynomial is constructed in terms of the values of the 
variable at the collocation points with the help of truncated Chebyshev 
series. Next, the first and second derivative of the variable are explicit- 
ly determined using the above interpolant such that (as an example we pre- 
sent only variable F(s) since the extension to other variables is very 
obvious) 
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HP N 

4L = Z A . F 

dc j k=0 Jk k 

(40) 

d 2 F N 

— = E B ik F k 
d£ j k=0 J K 

(41) 


j - 0,1, ... ,N 

where and are the elements of the derivative matrices and are given 
by Gottlieb et al . [1983] as 


with 


C 1 (-l) k+J 

A jk =■ J. iiJ — (j * k) 

r ^ ' 5k 

L k 


A. • = - v 

Jj 2(1 - c?) 


A _ 2N 2 + 1 _ _ a 

oo 7 NN 


C Q = C N , Cj - 1, (1 < j < N-l), 


(42) 


and 


B = AA = A 2 (43) 

It is clear that any higher derivative can be obtained by utilizing relation 
(43). 
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From our transformation, we obtain 


1 u 2 1_ = 4 

r 1 - 5* r 2 (1-s) 2 


and 


d _ dg d __2^ $ d 

dr dr d£ d£ 1 d£ 

d 2 d ,d, . d 2 c d 2 

dr 2 dr dr dc 2 


Here Si and S 2 are the scaling factors involved from the transformation. 

Writing the governing equations in terms of the new variable and evalu 
ating at the collocation points, we have 


Continuity: 


>1 . £ 


A jk F k + (3^77) F j + (3777) 6 J 


2 ' F, + fj!L ' 1 G, + = 0 


k=0 


(44) 


r-momentum: 


N 


3 2 i B-.F.+ f 

k=0 JK K l 1-e . 


Re u j] Si k i 0 A jk F k + i 1Re 


0) 


ReS x ^ 
d5 


i2Re n V . 


" J r 4(n 2 +l)-ji F r4(2n) 

i Tm- W ] j Wf + 


2( 12Re V.) N 

L] G + iReSl j A jk P k - a iRe W, F, 

1-e,- J k=0 JK K J J 


a 2 F . = 0 
J 


(45) 
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e -momentum: 


-C iReSj 


dV 

dC 


4 (2n) 2 < lRe v i> N 

+ U * i-] F, + S 2 I 6 1k G k 

j (l-«j) 2 J k=0 Jk 


N 2( iRen V .) 

+ \— — - Re U,] Si z A.. G k + f iReco - J 


l-C . J k=0 JK K 

J 


w j 


— eU -^ - 4 ^ n . 2 .t 1 .l] Si - 2 ( .--- n * Pi - oiRe U, G,- a 2 G. = 0 (46) 

1-5 j (l-5j) 2 J 1-?, J J J J 


Z-momentum 


-iRe Si 


F. + S 2 E B^ k H k + [_ RellJ 


k=0 




N 2(iRe n V.) 

s i 2 A -k. H. + [ iReai - 

k=0 JK K 3Z 


H- Rei“ 


4(n 2 ) 

d-^j) 2 


■] H. - iaRe W. H, - iaRe P. - a 2 H. = 0 

J J J J J 


(47) 


with the boundary conditions 


at 5 = -1 


F(-l) - G(-l) = H(-l) = 0 


for all n 


3P 

H 


= Normal momentum (48) 

5 = -1 evaluated at £ = -1 


at £ = 1 


If n = 0 


F(l) = G(l) = 0 

H(1 ) & P(l) are finite 


(49) 
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If n = ±1 


If M > 1 


F(l) ± 6(1) = 0 

H( 1) = P(l) = 0 (50) 

F(l) = 6(1) = 0 

H(1 ) = P(l) = 0. (51) 
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appendix a 


The three momentum equations plus continuity were discretized using a 
second order accurate central differencing. These equations are: 

Continuity 


(~) 0 \i+r F i-l) + — F i + — G. + «H. = 0 
2a r J 1 J 1 r J r J J 

j 


(1A) 


r j 


r-momentum 


C— ) 5.J * F J- 1 + F j + i ~ f j'-i . ("in + c.2) p .] 

° - 2 J J 


Re Ar 2 


2r . Ar r. 

J J 


1 [U. (J!I ti) + u' F,] + [« - n (J.) - . U ]F 

J 2Ar J J - J J 


r j 


2 [JJL. - (!i)] S, + Pjtl ' ?i - 1 = 0 , 

R e r j 2 r j 24r 


(2A) 


0 -moment urn 


(d) [ 6 W 26 J * G J-i t G “’ - G - 

Re Ar 2 2r . Ar 

J 


£1 ^±1 - + « 2 ) fi.] 

Oy Ay' v* 2 J 
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(3A) 


+ [i v : 


z- momentum 



, H . . - 2H . + H. . H. x1 - H. . 2 

pL) [ Jli J JiL + Ji I . (if_ + „s) h,j 

Re Ar 2 2r • Ar r . 2 ^ 

J J 

+ U, (Jll ill) + [- 1u + in (J.) + (i“) i + io W.] H. 

J 2Ar r . 3 r ^ J J 

j 

3W 

+ fi f — )j] Fj + iaPj = 0. (4A) 

3r 


The resulting matrices (or system of equations) are 4N x 4N, where N is the 
number of grid points involved. It is obvious that even for moderate values 
of N the memory requirement becomes very large. 

Once the code was operational, we found difficulty in obtaining a con- 
verged solution. Results of a sample calculation are given in Table 5. 

These results show that despite an initial decrease in the residual, it 
grows very fast and diverges quite rapidly. At first it was believed that 
the small numoer of grid points might have caused the problem out a suDstan- 
tial increase in the number of grid points showed that this was not the 
case. Next, attention was given to the accuracy of the method as a possible 
cause of the divergence. It was realized that the number of unknown factors 
involved (such as the number of grids, second order accuracy, boundary 
conditions, etc.) were too many to allow any direct checking of each factor 
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Table Al. Sample calculation showing convergence problem (M is the number of 
iterations.) 


M Residual 

1 51.03473794034 

2 15.83350551921 

3 11.31391553601 

4 8.397664066553 

5 5.701436944138 

6 4.314902731136 

7 2.065426550367 

8 1.209095939138 

9 .8214686963228 

10 .7912886143684 

11 .6251527839796 

12 .5471210004286 

13 .3662554796695 

14 .2757823153175 

15 .1838733412282 

16 .1157389438433 

17 .06495193935236 

18 .05683261464769 

19 .0617858871914 

20 .1043274836963 

21 .3839152251866 

22 .5645749241422 

23 1.967396398417 

24 6.97557243855 

25 10.17577382309 

26 21.72966773107 

27 59.67271919984 

28 151.0126448191 

29 437.9075245724 

30 924.9573854378 

31 1850.30198365 

32 4107.963293855 

33 12777.07307647 

34 40884.39728102 

35 71145.88056394 

36 119477.8852432 


separately. Nevertheless, we knew that it was essential to pinpoint the 
exact cause of the difficulty. Therefore, a model problem was chosen so 
that a quick analysis of some of the above unknowns was possible in a short 
time; and at the same time the problem had been solved by other workers so 
that comparison could be made. 

The temporal stability of Poiseuille flow in annulus was thought to be 
an ideal case for the above purpose. 
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APPENDIX B 


AXIAL VELOCITY PROFILE IN CONCENTRIC CYLINDERS 
The axial velocity for the Poiseuille flow in an annulus is given by 


W* = (Pi - Pa) [a 2 - r* 2 ♦ 0> 2 - a 2 ) t „ 
4uL £n(b/a) a 


(IB) 


where Pi, P 2 are the pressure, L is the length of the pipe, u is the dynamic 
viscosity, a and b are the inner and outer radius, respectively and r* is 
the dimensional radial distance. 

Define a new variable such that 


where 

and 

Therefore 


a r* = r(K-l) 
~ ~ 

K = — 
a 

2 


r 


2r* 

(b-a) 


then the velocity profile is transformed to 


or 


W* = ( p i - p 2_)_a 2 [1 _ 52 + K -.1, £ n g] 

4uL tn(K) 


U = — = [1 - K 2 + K2 - in 5]. 


Wo 


*n (K) 


(2B) 


(3B) 


(4B) 


(5B) 


41 


Furthermore, it would be advantageous if we normalized the above profile 
with respect to the maximum velocity occurring in the annulus. Therefore 


dW 

d? 


= - 2? + 


K 2 -l 
tn K 


1 

C 


the maximum velocity occurs where 


with 


— = 0 = - 2c m + 1 

dC inK C M 


S 2 


M 


K 2 -l 
tn K 2 * 


(6B) 


( 7 B) 


Then evaluating the velocity 


or 


which results in 


V«M> * 1 ’ C M * 777 *" 5 


M 


= 1 * + g £n C M 


(8B) 


W M^M> = 1 - 




+ K 2 - 1 

tn K 2 


tn £ 


2 

M* 


With more simplification we deduce 

= 1 + £n 
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Now normalizing the velocity with the above maximum value, we get 


— ■ W(5) - 
W M 


1 - C 2 + Cj, *n £ 2 
1 " + 5 M ln 


(10B) 


Since we want the non-dimensionalization to be with respect to half of the 
gap distance, that is 


r<_ii 


(b-a) (b-a) 


( 1 IB ) 


Then form (2B) we obtain the limits on c as 


1 < 5 < K. 


(12B) 
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